{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Crater attributes\n",
    "Prototype analysis for a single detected crater. Extract crater attributes like best-fitting ellipse, eccentricity, gradient, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "from pathlib import Path\n",
    "\n",
    "import numpy as np\n",
    "import rasterio\n",
    "from skimage.measure import EllipseModel, regionprops, label, find_contours, approximate_polygon\n",
    "\n",
    "# Get path to test images\n",
    "images_dir = Path.cwd().parent.parent / 'test/fixtures'\n",
    "\n",
    "# Construct filenames to grayscale image and individual crater annotations\n",
    "img_ctx_fpath = images_dir / 'img1.png'\n",
    "img_annot_fpaths = [images_dir / f'img1_annot{num}.png' for num in range(3)]\n",
    "grad_image_fpath = images_dir / 'diagonal_gradient.png'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/wronk/.virtualenvs/mars/lib/python3.7/site-packages/rasterio/__init__.py:219: NotGeoreferencedWarning: Dataset has no geotransform set. The identity matrix may be returned.\n",
      "  s = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)\n"
     ]
    }
   ],
   "source": [
    "orig_image = rasterio.open(img_ctx_fpath).read().squeeze()\n",
    "annot_images = [rasterio.open(img_annot_fpath).read().squeeze() \n",
    "                for img_annot_fpath in img_annot_fpaths]\n",
    "grad_image = rasterio.open(img_ctx_fpath).read().squeeze()\n",
    "\n",
    "label_objects = [label(annot_image) for annot_image in annot_images]\n",
    "props = [regionprops(label_object) for label_object in label_objects]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_crater_props(binary_label_img):\n",
    "    \n",
    "    if binary_label_img.ndim != 2:\n",
    "        raise RuntimeError('Image must be 2d for calculating region properties.')\n",
    "    \n",
    "    props = regionprops(binary_label_img)\n",
    "\n",
    "    print(f'Area = {props[0][\"area\"]}')\n",
    "    print(f'Major Axis = {props[0][\"major_axis_length\"]}')\n",
    "    print(f'Minor Axis = {props[0][\"minor_axis_length\"]}')\n",
    "    print(f'Eccentricity = {props[0][\"eccentricity\"]}')\n",
    "    print(f'Orientation = {props[0][\"orientation\"]}')\n",
    "\n",
    "    # Spatial moments?\n",
    "    # Euler number?\n",
    "    \n",
    "    # TODO: If we use EllipseModel, need to grab only the edge pixels (and not use filled blob)\n",
    "    coords = np.stack(np.nonzero(binary_label_img), axis=1)\n",
    "    em = EllipseModel()\n",
    "    em.estimate(coords)\n",
    "    print(f'Ellipse model params: {em.params}')\n",
    "    \n",
    "    \n",
    "def calculate_grad(image, good_inds):\n",
    "    \"\"\"Calculate the X and Y gradient for an image using numpy's gradient.\n",
    "    \n",
    "    Parameters\n",
    "    ==========\n",
    "    image: array-like\n",
    "        Image pixels.\n",
    "    good_inds: array-like\n",
    "        Bool mask for pixels to include when returninng mean gradient. \n",
    "        For example, pass the binary crater mask to only calculate the \n",
    "        gradient for pixels within the crater.\n",
    "    \n",
    "    Returns\n",
    "    =======\n",
    "    mean_h: float\n",
    "        Mean horizontal gradient of pixels specified by `good_inds`. Positive\n",
    "        is to the right.\n",
    "    mean_v: float\n",
    "        Mean vertical gradient of pixels specified by `good_inds`. Positive is\n",
    "        downward.\n",
    "    \"\"\"\n",
    "    # TODO: could improve function to take a list of good_inds and avoid\n",
    "    # recomputing gradient repeatedly\n",
    "    \n",
    "    if good_inds.dtype != bool:\n",
    "        raise ValueError('`good_inds` must be of type bool')\n",
    "        \n",
    "    np_mask = np.invert(good_inds) # Pixels we want to exclude should be True\n",
    "    \n",
    "    mean_h = np.ma.masked_array(np.gradient(image, axis=0),\n",
    "                                mask=np_mask).mean()\n",
    "    mean_v = np.ma.masked_array(np.gradient(image, axis=1),\n",
    "                                mask=np_mask).mean()\n",
    "\n",
    "    print(f'Grad horizontal:{mean_h}, Grad vertical:{mean_v}')\n",
    "    return mean_h, mean_v\n",
    "\n",
    "\n",
    "def calculate_wkt_border(binary_label_img, contour_tol=2.5, pad=1):\n",
    "    \"\"\"Pull out the polygon representation of the crater border.\"\"\"\n",
    "    \n",
    "    # TODO: apply and remove padding properly; add test\n",
    "    print(f'orig shape: {binary_label_img.shape}')\n",
    "    padded_image = np.pad(binary_label_img, pad_width=1, mode='constant')\n",
    "    print(f'padded shape: {padded_image.shape}')\n",
    "\n",
    "    contours = find_contours(padded_image, 0.5)\n",
    "    if len(contours) > 1:\n",
    "        warnings.warn(f'Found {len(contours)} contour objects. Expected 1.')\n",
    "    \n",
    "    coords = approximate_polygon(contours[0], tolerance=contour_tol)\n",
    "    closed = coords[0] == coords[-1]\n",
    "    \n",
    "    return coords, closed\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Object 0\n",
      "Area = 365\n",
      "Major Axis = 21.55794993684353\n",
      "Minor Axis = 21.55794993684353\n",
      "Eccentricity = 0.0\n",
      "Orientation = 0.7853981633974483\n",
      "Ellipse model params: [236.99999999996908, 44.00000000000697, 7.621885183733106, 7.621887402948319, 1.6470532156733175]\n",
      "Grad horizontal:0.8780821917808219, Grad vertical:-0.33013698630136984\n",
      "\n",
      "Object 1\n",
      "Area = 8166\n",
      "Major Axis = 143.50196551711082\n",
      "Minor Axis = 76.42045121336285\n",
      "Eccentricity = 0.8464052120638251\n",
      "Orientation = 0.0\n",
      "Ellipse model params: [174.4999999980865, 225.55693613864304, 27.124554096602242, 50.587078011778424, 1.5707963269431249]\n",
      "Grad horizontal:0.37447954935096744, Grad vertical:0.6970977222630419\n",
      "\n",
      "Object 2\n",
      "Area = 1129\n",
      "Major Axis = 50.04730585494444\n",
      "Minor Axis = 28.76355266027332\n",
      "Eccentricity = 0.818344999167847\n",
      "Orientation = -0.514135377558573\n",
      "Ellipse model params: [47.15534054087704, 21.710982033865484, 10.175457518522418, 17.68397594416444, 1.0568183030693015]\n",
      "Grad horizontal:-0.9694419840566874, Grad vertical:-1.8556244464127547\n"
     ]
    }
   ],
   "source": [
    "# Get some values from first crater\n",
    "for li, binary_label_img in enumerate(label_objects):\n",
    "    print (f'\\nObject {li}')\n",
    "    calculate_crater_props(binary_label_img)\n",
    "    calculate_grad(orig_image, binary_label_img.astype(np.bool))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "orig shape: (256, 256)\n",
      "padded shape: (258, 258)\n",
      "orig shape: (256, 256)\n",
      "padded shape: (258, 258)\n",
      "[[246.99803922 256.        ]\n",
      " [243.99803922 232.        ]\n",
      " [229.99803922 208.        ]\n",
      " [210.         192.00196078]\n",
      " [186.         184.00196078]\n",
      " [160.         185.00196078]\n",
      " [133.         197.00196078]\n",
      " [113.00196078 219.        ]\n",
      " [104.00196078 244.        ]\n",
      " [105.         256.99803922]\n",
      " [246.99803922 256.        ]]\n",
      "orig shape: (256, 256)\n",
      "padded shape: (258, 258)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3EAAAEZCAYAAAA5TKxSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de5RcVZmw8edNJyEJ15BAuCQhAaKIIrfIRRzAoCMwKKCSARER0DiKozCyFHEhYcQZxlFY6Ix+hpEBB0dAnUhAAQVUFAVJMNyJRAFJDIRLJMg9nf39UafSlU53ujtdp0+dOs9vrVpdZ/fpOm8qye791t7n3ZFSQpIkSZJUDsOKDkCSJEmS1H8mcZIkSZJUIiZxkiRJklQiJnGSJEmSVCImcZIkSZJUIiZxkiRJklQiuSVxEXFoRCyKiMURcWZe15EkSSo7x02SBiLy2CcuIjqA3wNvB5YAdwDHpZTub/rFJEmSSsxxk6SBymsmbh9gcUrpjymlV4ArgCNzupYkSVKZOW6SNCB5JXHbA481HC/J2iRJkrQ2x02SBmR4UReOiFnArOxw76LikJSflFIUHYMktQvHTlL76+/YKa8kbikwqeF4Yta2RkppDjAHICKaf2OeJElSOfQ5bgLHTpK65LWc8g5gWkRMjYiRwLHAvJyuJUmSVGaOmyQNSC5JXEppFfBx4AbgAeCqlNJ9eVxLkgbCMt6SWo3jJkkDlcsWAwMOwiUBUltqtXviLOMtqV04dpLaU3/HTrlt9i1JLcgy3pIkqfQKq04pSQXoqYz3vuv7AT/tllrf3nsPrFDjggULnkopbZVTOJKUO5M4SeqmWxlvSS1u/vz5Azo/Ih7NKRRJGhImcZKqxDLekiSp9LwnTlKVWMZbkiSVnjNxkiojpbQqIuplvDuASyzjLUmSysYtBiTlptW2GNgQ9k9S6xvoWCYiFqSUpucUzpCwb5Lak1sMSJIkSVIbMomTJEmSpBIxiZMkSZKkEjGJkyRJkqQSMYmTJEmSpBIxiZMkSZKkEjGJkyRJkqQSMYmTJEmSpBIxiZMkSZKkEjGJkyRJkqQSGV50AJIkSWptmwETs0cH8BOgs9CIpGoziZMkSRJjgKOAacCZxx/PiMcfJ5YurT2ee26tc6/Pzn156MOUBERKqegYiIjig5DUdCmlKDqGwbJ/klrfQMcyEbEgpTQ9p3CGRLP7ph2BB6ZOZeTDD/f7Z5474AC2uvVWEzmpifo7dvKeOEmSpAp7M7Bo7Nj1J3CjRsG0abDPPmuaNr31Vpbvvz8b5R+ipG5M4iRJkirqvrPP5taNNmL4ihW1ho02gjPOgG98A665BhYuhKeeghdegN//Hm6/HWbPXvPzm/3mNyzbd18TOWmIuZxSUm5cTilpKLiccuCGRXDf8cezy+WXdzVutRVcfTXsv3/fL3DuuWslc0/tvTeT77yTF1tgXCmVmcspJUmStI7RHR0sPvDAtRO4XXaB227rXwIHcM45ayVx4xcs4OHdd3dGThoiJnGSJEkVsSmweOedmfqLX3Q1zpgBv/417LjjwF6sWyI3YeFCfggmctIQMImTJEmqiIuGD2e7RYu6Gk46Ca67DsaO3bAX7JbIHQomctIQMImTJEmqgD023pgTV6/uajjvPPjWt2DkyMG9cC+JnJsRS/kxiZMkSaqAb0+dyrB6EjdjBnzucxBNqj/VQyL3nua8sqQemMRJkiS1uYPGjmW3e+/tavjiF5t/kXPO4fFddmn+60pah0mcJElSG9t66625ZOLEroZ3vhP22y+Xa4184YU1z/+YyxUkgUmcJElSWzttv/3Y8Z57uhrOOy+fC6XEJsuXrzk0iZPyYxInSZLUpiZOnMgJjdUojz0W3vjGfC729NOMfOklAFYCT+dzFUmYxEmSJLWliRMn8p/veQ8T60lcRwece25u13v81lvXPHcWTsqXSZwkSVIb2mP33XnX7bd3NZx0ErzmNbld708///ma5yZxUr5M4iRJktrM5MmTOX3aNLjttlrDyJFw9tm5XnPjJ55Y89wkTsqXSZwkSVKb2XHKFGbcfHNXw0c/CpMn53pNkzhp6JjESZIktZHJkyfzb9Onw9131xrGjIHPfjb365rESUPHJE6SJKmNbLvVVuxzzTVdDaedBhMm5HrNe+65h87Fi9cc/yHXq0kyiZMkSWoTEydO5NK3vhUeeqjWsPnmcMYZuV932aOPsvXLLwOwGng09ytK1TZ8MD8cEY8AzwGdwKqU0vSI2BK4EpgCPALMTCmtGFyYkiRJ5Zf32GncJpuwy5VXdjV8+tMwduzggu6H0cuXr5kZeAx4NfcrStXWjJm4t6aU9kgpTc+OzwRuSilNA27KjiVJklSTy9hpm2224Ufvehc89litYeut4ROfaEa8fRq9bNma594PJ+Uvj+WURwKXZc8vA47K4RqSJEntoiljpy1HjmT7Sy/tajjrLNhkk8HG1qe77rqLH1144Zpjkzgpf4NN4hLwk4hYEBGzsrYJKaX6xzGPA/neSStJklQeuY2dPrhyJSxfXjuYNAk+8pHBxtovL6xYwfFPP73m+P4huapUbYO6Jw54S0ppaURsDfw0Ih5s/GZKKUVE6ukHs45rVk/fk6TB8p5dSS0ql7HTzuPG8alVq7oaPv95GDWqeVH34v777+eGI45gdna8Arg096tKGtRMXEppafZ1OTAX2Ad4IiK2Bci+Lu/lZ+eklKY3rAeXpGbznl1JLSWvsdOpL73EsGefrR1MmwYnnpjPH6DBH/7wB47cf39Of/75NW2zgWdyv7KkDU7iImLjiNi0/hz4W+BeYB5Q7zlOBK4ebJCS1CTesyupMHmNnbYGPvzCC10N//zPMGJEEyJev87OTj61ciWbZ8cPAl/P/aqSYHDLKScAcyOi/jr/m1K6PiLuAK6KiFOobRMyc/BhStKA1e87ScA3U0pz8J5dScXKZew0e+RINn7lldrBG98IM4dm6DXywQf5cMPxPwGrejtZUlNFSj0uux7aIHpZ+y2p3FJKUdS1I2L7xvtOgH8E5qWUtmg4Z0VKaZ0NlLrdd7L3kAQsaYMNdCwTEQvKfjtHfew0EngWWHP327x58M535n79JY89xkOTJ/PW7Ph64LDcryq1v/6OnfLYYkCSCuc9u5KqYBMaErjhw+Hww4fkuqN+8pM1CdwqarNwkoaOSZyktuM9u5Kq4hng8cg+uF+1ChYvzv2aTy1dyl8+9KE1x98AHsj9qpIaDXaLAUlqRd6zK6ky7u3oYJv69gJ33w2vfW2u1xv5zW+yc/Z8BazZXkDS0DGJk9R2Ukp/BHbvof1p4JChj0iS8nN/Rwdvqydxd90FxxyT27WW33svG33hC2uOZ+OWAlIRXE4pSZJUYre//HLXwd1353adZcuW8X+77eaWAlILcCZOkiSpxNZK2+66K5drdHZ28o7ttuN3DW1uKSAVx5k4SZKkEnsQWDMX96c/wV/+0vRrPHLvvXwN6MiOrweua/pVJPWXSZwkSVKJrQIWdXR0NTRxSeW9997Lw//938Qee3BQw/XcUkAqlkmcJElSyf2us7ProElJ3J233MLPd9uNqSefzI4N7efjlgJS0UziJEmSSm6tO+HuvHPwr3fRRWxx0EF8vKFtBXACcPagX13SYJnESZIklVxjEpcuuwyuvHKDXuemq6/m4b/7O3Y/7bS1Zt/mAbsClw8iRknNY3VKSZKkkvsFsBDYA4jVq0nvex8B8Pd/3+/XuPWLX2Sns89mSkpr2lYAn8DkTWo1JnGSJEkl1wm8A7gZeD39T+SuvfZalv/xj+w7dy4H/Pzna31vHvAR4PGcYpa04SI1fNpSWBARxQchqelSSlF0DINl/yS1voGOZSJiQUppek7hDIne+qat6UrkAFZHcMPMmTy20049vs5DF1/M+U8+SUNtS2ffpAL1d+xkEicpNyZxkoaCSdzauidyA+Hsm1Ss/o6dLGwiSZLURpYDM4D7BvhzJwBHYgInlYH3xEmSJLWZeiL3ZeDAXs7ZoeH53wNX5R2UpKYxiZMkSWpDy4EP9HHO/sBvhiAWSc3lckpJkqSKMoGTyskkTpIkSZJKxCROkiRJkkrEJE6SJEmSSsQkTpIkSZJKxCROkiRJkkrEJE6SJEmSSsQkTpIkSZJKxCROkiRJkkrEJE6SJEmSSsQkTpIkSZJKxCROkiRJkkrEJE6SJEmSSsQkTpIkSZJKxCROkiRJkkrEJE6SJEmSSsQkTpIkSZJKxCROkiRJkkrEJE6SJEmSSsQkTpIkSZJKxCROkiRJkkqkzyQuIi6JiOURcW9D25YR8dOIeCj7OjZrj4j4akQsjoi7I2KvPIOXJElqNY6dJOWtPzNxlwKHdms7E7gppTQNuCk7BjgMmJY9ZgHfaE6YkiRJpXEpjp0k5ajPJC6ldAvwTLfmI4HLsueXAUc1tH871dwGbBER2zYrWElq5KfdklqRYydJedvQe+ImpJSWZc8fByZkz7cHHms4b0nWJkl5uBQ/7ZZUDo6dJDXNoAubpJQSkAb6cxExKyLmR8T8wcYgqZr8tFtSGTl2kjRYG5rEPVEf/GRfl2ftS4FJDedNzNrWkVKak1KanlKavoExSFJP/LRbUity7CSpaTY0iZsHnJg9PxG4uqH9A9m9J/sBzzYMpiRpSPlpt6QW4thJUtMM7+uEiPgucDAwPiKWAOcA5wNXRcQpwKPAzOz0HwOHA4uBF4CTcohZktbniYjYNqW0bDCfdgNzACJiwEmgpA23MbX/nH8BngZWFRvOBnHsJClvUfuguuAgHCRJbSmlFHlfIyKmANemlN6QHf878HRK6fyIOBPYMqX06Yj4O+Dj1AZL+wJfTSnt04/Xt3+SchQRjIzgHcD7UuJdKTG64ftps81g/HjSuHEwfjxsuSVp/HgYNw7GjSONG8fwffeFKVMGcs0FZV+SaN8ktaf+jp1M4iTlJu8krvHTbuAJap92/xC4CphM9ml3SumZiAjgP6hVs3wBOCml1OdySfsnqfk233xzSIm9Ozs5b5dd2P/RR+Gppzb8BTs64IYb4JBD+nW6SZykVmUSJ6lwQzETlzf7J6k5dthhByKC12+0EdcedxxcfjksXtzzyZMmwYsvwjPPwOrV/bvA6afDBRf061STOEmtqr9jpz7viZMkSdpQu+66K5NGj2bucccx+vvfh9tug9mz1z1x4kQ4/vjaY7fdam2rV8Nf/lKbpXv66drXxue/+EXt9SSpYkziJElSU02ZMoU377kne//5z3xo9Gg2+9WvYMGCdU/cbDN473vhhBPgwANhWLei2cOGwZZb1h49+fKXu5K4KP3EvyT1m0mcJElqiq3Hj+ecgw/msKeeYuqNN8Jzz6170vDhcPjh8P73wxFHwOjR657TX41LLbsngJLUxkziJEnSoOw5fDhzDjyQ1y1cyMbf/37PJ735zbXE7ZhjalUmm8EkTlJFmcRJkqQB2w54H3DW5MmM/dOf4Oab1z1p2rRa4nb88bDTTs0PwiROUkWZxEmSpH4bDnwV+IcIIiX405/WPmGrreDYY2vJ25velO+9ao0Vtk3iJFWISZwkSeqXMcD3gMNh7QRq1Cg46qhagZK3vx1GjBiagJyJk1RRJnGSJKlPWwILtt2WKcuWdTUedBCcdBIcfXSt0uRQa0zirE4pqUJM4iRJ0nrtNGIEt4wZw3aNCdzZZ8O55xabPDkTJ6miTOIkSVKPIoKZb3gDV6xYAUuW1Bvha1+DU08tNjgwiZNUWSZxkiSpR6fsuisXL1kCK1bUGkaMgMsvh5kziw2szsImkirKJE6SJK3jrD324IuLFsGLL9YaNtkEfvhDOOSQYgNr5EycpIoyiZMkSWscffTRHLZ8OR+67Tbo7Kw1brUVXHcd7L13scF1Z2ETSRXlx1aSJAmAD37wg1y+++58+NZbiXoCN3Uq3Hpr6yVw4EycpMqyx5MkSfzDrFn8x6hRjJk9u6tx991rCdy0aYXFtV4mcZIqyuWUkiRV3EdOPpkLnn6a0T/4QVfjQQfB1VfD5psXF1hfLGwiqaJM4iRJqrAPHnMMX330UUbedFNX47vfDd/5DowaVVxg/eFMnKSKMomTJKmiZh5yCP/18MN0zJ/f1ThrFnz969DRUVxg/WVhE0kVZRInSVIFTQK+cvvtdPz1r12NZ58N555bnoTImThJFWUSJ0lSxbweuHnkSLauJ3AR8LWvwamnFhrXgJnESaookzhJkirkzcCPOzrY/JVXag0jRsDll8PMmYXGtUEsbCKpokziJEmqiMOBHwwbxqj6HnCbbAI//CEcckihcW0wZ+IkVZQ9niRJFfABYF4Eo+qJz1Zbwc9/Xt4EDixsIqmyTOIkSWpz/whcBnTUlx9OnVrbxHvvvYsMa/CciZNUUfZ4kiS1sR2ACxobdt+9lsBNm1ZQRE1kEieporwnTpKkNvYJGn7Z77UX3HwzbL55gRE1kYVNJFWUPZ4kSW1qU+BDjQ3nndc+CRw4EyepsuzxJElqU6cAm9UPXvc6eMc7CowmBxY2kVRRJnGSJLWhDuCTjQ2nndZ+s1XOxEmqKHs8SZLa0HsimFI/GDcOTjihwGhyYhInqaLs8SRJakNf2X77roOPfQxGjy4umLxY2ERSRdnjSZLUZg4ZM4aJS5bUDkaOrCVx7ciZOEkVZY8nSVKb+dcJE7oO3vc+2Gab4oLJk4VNJFWUSZwkSW1k+vjx7P3II10Np59eWCy5cyZOUkXZ40mS1EbO3HhjhtXvFXvb2+CNbyw2oDyZxEmqKHs8SZLaxKbAofV74QD+6Z8Ki2VIWNhEUkXZ40kqrYi4JCKWR8S9DW2zI2JpRCzMHoc3fO+zEbE4IhZFRJvteizBGVtswcadnbWDdtzcuztn4iRVVJ89noMkSS3sUuDQHtovTCntkT1+DBARuwLHAq/PfubrEdExZJFKOesATnn++a6Gdtzcu7sWLWzi2ElS3vrTu1+KgyRJLSildAvwTD9PPxK4IqX0ckrpYWAxsE9uwUlD7MRNN2X7V1+tHbTr5t7dte5M3KU4dpKUoz57PAdJkkro4xFxd/Zp+NisbXvgsYZzlmRtUlv46EsvdR206+be3bVoEufYSVLeBtPjOUiS1Iq+AewE7AEsA74y0BeIiFkRMT8i5jc7OCkP+wHT67Nw7by5d3flK2zi2ElSU2xoj+cgSVJLSik9kVLqTCmtBi6m6xPtpcCkhlMnZm09vcaclNL0lNL0fKOVmuOMxgSmnTf37q5FZ+J64dhJUtNsUI/nIElSq4qIbRsOjwbqhQXmAcdGxEYRMRWYBvx2qOOTmm0H4KjGZKadN/furkULm/TEsZOkZtqgJM5BkqRWEBHfBX4DvDYilkTEKcCXIuKeiLgbeCtwOkBK6T7gKuB+4Hrg1JRSZ0GhS03zCWqVKYH239y7uxLNxDl2ktRMw/s6IRskHQyMj4glwDnAwRGxB5CAR4CPQG2QFBH1QdIqHCRJylFK6bgemr+1nvO/CHwxv4ikobUp8KHGhnbf3Lu7Fk3iHDtJylukxpuCiwoiovggJDVdSqm11zf1g/2TWtlpwIX1g112gfvua6lkJnd/8zfwq1/Vnt9yS+24HyJiQdmXJNo3Se2pv2OnCvX0kiS1jw7gk40Np59erQQOWnYmTpLyZo8nSVIJHQ1MqR9UZXPv7kpU2ESSmskkTpKkElprFq4qm3t350ycpIqyx5MkqYT2bzyoyube3ZnESaooezxJkkpmFF3bCqRRo6qzuXd3jcXZTOIkVYg9niRJJTOm8aCKyyjrnImTVFH2eJIklUxjEhdjxvR6XtsziZNUUfZ4kiSVzFppm0lcjdUpJVWISZwkSSWz1gJKk7gaZ+IkVYg9niRJJeNMXMbCJpIqyh5PkqSSMYnLOBMnqaLs8SRJKhmTuIxJnKSKsseTJKlk3GIgY2ETSRVlEidJUsk4E5dxJk5SRdnjSZJUMiZxGQubSKooezxJkkrGJC7jTJykirLHkySpZEziMiZxkirKHk+SpJJZq5SJhU1qLGwiqUJM4iRJKhln4jLOxEmqKHs8SZJKxiQuY2ETSRVljydJUsmYxGWciZNUUfZ4kiSVjElcxiROUkXZ40mSVDImcRkLm0iqKJM4SZJKZq20zeqUNc7ESaoQezxJkkrGmbiMhU0kVZQ9niRJJWMSl3EmTlJF2eNJklQyjQson3/11cLiKJxJnKSKsseTJKlkljU87/jtbwuLo3AWNpFUUSZxkiSVzA8bng+fO7ewOArnTJykirLHkySpZL7X8Lzjxhth5crCYimUhU0kVZQ9niRJJfMH4M7sebzyCsybV2Q4xXEmTlJF2eNJklRCjbNx6Xvf6/W8tmYSJ6mi7PEkSSqhtdK266+v5pJKC5tIqiiTOEmSSqjySyob74cDkzhJlWISJ0lSSa01G1e1JZWNSVyESZykSjGJkySppNZK2264oVpLKr0fTlKF2etJKqWImBQRP4uI+yPivoj4ZNa+ZUT8NCIeyr6OzdojIr4aEYsj4u6I2KvYP4E0eI1LKnn5ZbjmmgKjGWImcZIqzF5PUlmtAj6VUtoV2A84NSJ2Bc4EbkopTQNuyo4BDgOmZY9ZwDeGPmSp+daajbvqqqLCGHoWNZFUYSZxkkoppbQspXRn9vw54AFge+BI4LLstMuAo7LnRwLfTjW3AVtExLZDHLbUdJVdUulMnKQKs9eTVHoRMQXYE7gdmJBSWpZ963FgQvZ8e+Cxhh9bkrVJpVbZJZWNhU1M4iRVTJ+9nvedSGplEbEJ8APgtJTSWlMQKaUEpB5/cP2vOSsi5kfE/CaFKeWqkksqW3gmzrGTpLz1p9fzvhNJLSkiRlBL4L6TUvq/rPmJ+jLJ7OvyrH0pMKnhxydmbetIKc1JKU1PKU3PJ3KpuRqTuM7rrqvGksoWTuJw7CQpZ332et53IqkVRUQA3wIeSCld0PCtecCJ2fMTgasb2j+QfeK9H/Bsw7JLqdQal1R2vPpqNZZUtnBhE8dOkvI2oI+umnnficuVJA3SAcAJwIyIWJg9DgfOB94eEQ8Bb8uOAX4M/BFYDFwMfKyAmKXcVG7j79aeiVvDsZOkPAzv74nd7zuJhk+9UkopIgZ030lKaQ4wJ3vtAd+zIqnaUkq/Anr7+P2QHs5PwKm5BiUV6HvAv2bPV/3oRwxfuRI226zIkPJVgsImjp0k5aVfvV5e951IkqTmaFxSOXzVqvZfUtniM3GOnSTlqT/VKb3vRJKkEmhcRPnQv/xLYXEMiRZO4hw7ScpbpLT+2fiIeAvwS+AeoN5jnkVtbfdVwGTgUWBmSumZrOP6D+BQ4AXgpJTSetduuyRAak8ppdaqNrAB7J9UJjtRu+kT4NWODjqeeophW2xRZEj5+fOfYfvstrFttoFl/c95ImJBntVnHTtJ2lD9HTv1mcQNBTsiqT2ZxElDbwFQ32TsdzNm8Ibrr2fEiBFFhpSPJUtgUrYCcbvtYGn/Vx/mncQNBfsmqT31d+zUWusPJEnSoPxXw/M9b76ZB2fPLiqUfJWgsIkk5cVeT5KkNvL/qN1gVbfLv/87L952W1Hh5KeF74mTpLzZ60mS1EYStQ0UH8yOR7z6Kh3vfS88/XSBUeXAJE5ShdnrSZLUZlYCR2VfAUYuXUrnMcfAqlUFRtVkJnGSKsxeT5KkNrQIeH/DccfPfsbqz3ymqHCarzGJi9LXUJKkATGJkySpTV0DzG44HnbBBXDFFQVF02QWNpFUYfZ6kiS1sX9m7UInnHwy3HVXQdE0kcspJVWYvZ4kSW2se6ETXnwRjjqq/IVOTOIkVZi9niRJba57oRMeeQSOPbbchU5M4iRVmL2eJEkV0L3QCTfeSOenP11QNE1gYRNJFWYSJ0lSRXQvdNJx4YXw3e8WFM0gWdhEUoXZ60mSVCHdC52sLmuhE5dTSqowez1Jkiqke6GTYS+9xKtHHFG+QicmcZIqzF5PkqSK6V7oZMSSJbxw5JHlKnRiEiepwuz1JEmqoO6FTsbceisrTz21qHAGzsImkirMJE6SpIrqXuhkszlzWPnNbxYUzQBZ2ERShdnrSZJUYd0LnYz82MeY9aY3sbpxpqsVuZxSUoXZ60mSVGHdC52MWr2ai+bPZ84228ALLxQYWR9M4iRVmL1eHwLYbj2PbYoLTZKkpuhe6GQ08A9PPsmfN98crrxy7aWLrcIkTlKF2eutx57AA8DS9TyWAXcCOxUUoyRJzbAIOBiY39C23apVcOyxpIMOYvWddxYTWG8sbCKpwkzierEncMvIkby2n+f+DBM5SVK5/Q7YBzgZeKKhPX75S5g+nc5Zs+Cpp4oJrjsLm0iqMHu9HkwfNoxbRoxgk1deqTVstBFst13PjxEjAJiEiZwkqfwS8N/Aa4CvAK9m7cNSouPii+ncaSdeueCC4veUczmlpAqz12sQEXxor724Y/PN2eTV7NfW2LHw61/D0qU9P66/HkaPBkzkJEntYyVwBrAbcH1De8fKlYz81Kfo3H13uPnmYoIDkzhJlWav12Dmzjtz8cMPw4oVtYaxY+HGG2GvvXr/oRkz4NprTeQkSW1pEXAY8E5gcUN7x/33wyGHkN79bnjkkaEPzCROUoXZ62X2jmDOQBO4OhM5SVKbuxZ4PfAZ4K8N7TF3LqumTYPPfx6ef37oArKwiaQKM4mjlsDd3NHBZvX1/QNJ4OpM5CRJbe4V4EvU7pf7dkP78FWr4AtfYOX228MVVwzNlgQWNpFUYZXv9QL4/rBhg0vg6npI5C5rWqSSJLWGZcCJwP6svSXBZs8+C8cdx7LXvAYWLsw3CJdTSqqwyvd6k4ApnZ21g1GjNjyBq8sSufqvlgOAjQYZoyRJreg2et6SYNvFi+ncc0/mTpjAwptuyufiJnGSKsxer9HWWw8ugaubMYNXBv8qkiS1vN62JOgAjl6+nClvfzuzx41jxZNPNvfCJnGSKsxeT5IkDVpvWxJskRKzn3mGJ7bZhotHjuS8LbaA730P7rgDnnxyw++fs7CJpAobXnQAkiSpfdS3JDgCuBDYOWvfZfVqdlm9Gp59FmbOXHP+y8OHM2LaNNhhB1L2GL7zzjB1KkyZAuPG9XUl+swAAAqZSURBVJykWdhEUoWZxEmSpKa7FvgJcBpwNrBJL+dttGoVPPBA7dGDNGYMnZMns3rSJDonTWL15Ml0TprE8IULGVM/ySROUsVEGooywH0FEVFYEJOBR9ccTIZHH13P2f33UgSjsuejgJeb8qpSuaSUclvjFBGTqFU5n0Dttpw5KaWLImI28GGgfgPOWSmlH2c/81ngFKAT+ERK6YZ+XKf4TlIquXHAgcCUbo+pwKbNusgAxjMRsSClNL1Zly6CfZPUnvo7dnImTlJZrQI+lVK6MyI2BRZExE+z712YUvpy48kRsStwLLX9ircDboyI16SUOoc0aqmCngbm9vK9sayb3NUfTU3yJKmNVD6JW+tjrMcfh1/9Ct7ylkG95neOPZbje7uGpKZIKS2jtl0VKaXnIuIBYPv1/MiRwBUppZeBhyNiMbXq6L/JPVhJvVqRPX7Xy/d7S/LelXdgktTCKr+I/DHg/vrBK6/AoYfWErkNdMm7380xV1655vhmcLsBKWcRMQXYE7g9a/p4RNwdEZdExNisbXtq/+XrlrD+pE9SC6gneHOpFUr5JLVPZILaNgZbAyxbVlR4klSIyidxAO8BHq8fPP/8hidy11zD++fOZWR2uBg4sSkRSupNRGwC/AA4LaW0EvgGsBOwB7WZuq9swGvOioj5ETG/qcFKaqrVZDe/brNNwZFI0tAyiQMeBN4KPFEvYbwhidw117DqqKPWSuDeSu2jfkn5iIgR1BK476SU/g8gpfRESqkzpbQauJjakkmApcCkhh+fmLWtI6U0J6U0veyFDyRJUnsyics8CByc0gYlche97W28euSRDM82HjWBk/IXEQF8C3ggpXRBQ/u2DacdDdybPZ8HHBsRG0XEVGAa8NuhileSJKlZ+kziImJSRPwsIu6PiPsi4pNZ++yIWBoRC7PH4Q0/89mIWBwRiyLiHXn+AZqpnsgNZGnlhTNm8NGbbmJEVtrYBE4aMgcAJwAzuvVDX4qIeyLibmr/HU8HSCndB1xF7TbY64FTrUwpKQ9VGjtJKkaf+8Rln2pv21jGGzgKmAn8tZcy3t+ltoRpO+BGYL1lvFttr5NdgJ8B9RX2ncDLPWwkOiabeaszgZPWluc+cUOl1fonSesa6J63ee8TV8Wxk6TmaNo+cVUs412/R66eyHWwbsLWExM4qS39FVhUdBAbYDzwVNFBDFAZY4Zyxl3GmKGXuCMG/HnRDk2JphdVHDtJGloD2ieuWxnvA6iV8f4AMJ/aprsrqHVStzX8WCnLeNcTuf8B+vNR3ccwgZPa1KIyFjiJiPlli7uMMUM54y5jzFDOuKs0dpI0dPpd2KTZZbzLUML7QeBNwCbreZwCvA6YU1CMkiSpNVVx7CRpaPRrJq63Mt4N378YuDY77FcZ75TSHLLcp9XXdT+/nu9dMmRRSJKksqj62ElSvvpTndIy3pKqrKwT7WWMu4wxQznjLmPMUJK4HTtJylt/qlO+BfglcA9Qr+5xFnActeUACXgE+Eh2Iy8R8TngZGAVtSUE1/VxDT9NktpQO1SnlKSBcuwkaUP1d+zUZxI3FOyIpPZkEidJ+XDsJLWnpm0xMESeonbrWRnLHdeVtVxznfEXqx3jz7WE91CIiEOBi6jtNPJfKaXzCw6pRxHxCPActW0tV6WUpkfElsCVwBRqn/jPzKrgFSYiLgGOAJanlN6QtfUYZ7Yc7SLgcOAF4IMppTtbJObZwIeBJ7PTzkop/Tj73mep1bzqBD6RUrqhgJgnAd8GJlCb8ZmTUrqoBO91b3HPpoXf7wI5diqe8RerHePv99ipJWbioJxlgxsZf7GMv1hlj78nEdEB/B54O7Vy33cAx6WU7i80sB5kSdz0lNJTDW1fAp5JKZ0fEWcCY1NKnykqxiymA6ntu/fthoSoxzgj4nDgH6klFvsCF6WU9m2RmGfTpA2bc4q5t42mP0hrv9e5b5Ddbsre9xp/sYy/WIONv99bDEhSxewDLE4p/TGl9ApwBbUNecviSOCy7Pll1AbDhUop3QI80625tziPpJY4pZTSbcAW3YpCDIleYu7Nmg2bU0oPA/UNm4dUSmlZfSYtpfQcUN9outXf697i7k1LvN+SVASTOEnq2fbAYw3Hrbz5bgJ+EhELImJW1jahXjABeJzaErVW1Fucrf7+fzwi7o6ISyJibNbWcjF322i6NO91t7ihJO+3JA2VVkriSlE2eD2Mv1jGX6yyx192b0kp7QUcBpyaLQFcI9XWzbfG2vn1KEucDHLD5qHSw0bTa7Tye93sDbLbXNn7XuMvlvEXa1Dxt0wSl21gWVrGXyzjL1bZ4+9FvzbfbQUppaXZ1+XAXGpLyp6oL4nLvi4vLsL16i3Oln3/U0pPpJQ6U0qrgYvpWsLXMjH3tNE0JXive9sgu9Xf76KUve81/mIZf7EGG3/LJHGS1GLuAKZFxNSIGAkcS21D3pYSERtnRSCIiI2Bv6W2gfA84MTstBOBq4uJsE+9xTkP+EDU7Ac827AUsFCtvmFzbxtN0+LvtRtkS9IApJQKfQCHAouo3ZB8ZtHx9DPmR6ht4LkQmJ+1bQn8FHgo+zq26Di7xXwJtU9d721o6zFmIICvZn8ndwN7tWj8s6l96rowexze8L3PZvEvAt5RcOyTgJ8B9wP3AZ8s0/u/nvhL8f4P8s9+OLUKlX8APld0PL3EuCNwV/a4rx4nMA64Kfv3dSOwZQvE+l1qy+FepXb/0im9xZn9P/jP7L2/h1r1zVaJ+X+ymO6mlkhs23D+57KYFwGHFRTzW6gtlby78f9nCd7r3uJu6fe7oPfKsVP+8TpuKjZ+x059XaPgP2BH1vnuCIykNgjZteh/OP2I+xFgfLe2L9U7UuBM4N+KjrNbfAcCe3X7z9xjzNkvzeuy/xD7Abe3aPyzgTN6OHfX7N/SRsDU7N9YR4Gxb1vvTIBNqSUFu5bl/V9P/KV4/3348OGjnR6OnYYsXsdNxcbv2KmPv4Oil1OWvYR3o5Yr590olbC0d6Ne4u9NS5WdTiUt9123nvh701LvvyS1GcdOQ8BxU7G/tx079f13UHQSV9bywGUu592oNOWm16NUZafLWu67zrLfklS4svax7TB2Kt3v7R6U7ve2Y6eeFZ3ElVVblPNuVMaYKVnZ6bKW+66z7LckaRDaauxUtngzpfu97dipd0UncaUsD5zKXc67UcuXm16fVKKy02Ut911n2W9Jahml7GPbZOxUmt/bPSnb723HTuuPv+gkrhQlvBu1QTnvRi1dbrovZSk7XdZy33WW/ZakluLYqTil+L3dmzL93nbs1PffwfDmhjwwKaVVEfFx4AZq1ZYuSSndV2RM/TABmFv7u2E48L8ppesj4g7gqog4BXgUmFlgjOuIiO8CBwPjI2IJcA5wPj3H/GNqVX4WAy8AJw15wN30Ev/BEbEHtan0R4CPAKSU7ouIq6iVdV0FnJpS6iwi7swBwAnAPRGxMGs7i/K8/73Ff1xJ3n9JahuOnYaG46bCf287durj7yCyspaSJEmSpBIoejmlJEmSJGkATOIkSZIkqURM4iRJkiSpREziJEmSJKlETOIkSZIkqURM4iRJkiSpREziJEmSJKlETOIkSZIkqUT+P4WBJK2Kx1acAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x504 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, axes = plt.subplots(ncols=3, figsize=(15, 7))\n",
    "plt.gray()\n",
    "\n",
    "for ii, ax in enumerate(axes):\n",
    "    ax.imshow(annot_images[ii])\n",
    "    contour, closed = calculate_wkt_border(annot_images[ii])\n",
    "    \n",
    "    ax.plot(contour[:, 1], contour[:, 0], '-r', linewidth=3)\n",
    "    if ii == 1:\n",
    "        print(contour)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
